function [source_recon_results ] = osl_inverse_batch( S )

% [source_recon_results ] = osl_inverse_batch( S )
%
% Run LCMV beamformer in PCA reduced sub space using SPM Beamformer toolbox
%
% Mark Woolrich 2012

global OSLDIR;   % added by DM

% Required settings:
try, S.D_continuous=S.D_continuous; catch, error('S.D_continuous not specified'); end; % SPM MEEG object file path of continuous data to run the analysis on
% AND/OR
try, S.D_epoched=S.D_epoched; catch, S.D_epoched=[]; end; % SPM MEEG object file path of epoched data to run the analysis on

try, S.time_range = S.time_range; catch, error('S.time_range not specified'); end; % Time range (from to) within trial, in secs, need to all be the same duration and one for each condition
try, S.dirname = S.dirname; catch, error('S.dirname not specified'); end; % dir name to store results in
try, S.mni_coord = S.mni_coord; catch, error('S.mni_coord not specified'); end; % num_voxels x 3 matrix of mni coords

% Optional settings:
try, S.epochinfo=S.epochinfo; catch, S.epochinfo=[]; end; % epochinfo for passing to spm_eeg_epoch (see help spm_eeg_epoch). if empty, will look inside D_epoched for D_epoched.epochinfo
try, S.freq_range = S.freq_range; catch, S.freq_range = []; end; % Frequency range (from to) in Hz
try, S.conditions=S.conditions; catch, S.conditions={'all'}; end; % list of triggers (conditions) to include in the analysis, e.g. trigger={'Motorbike','Neutral face'};
try, S.modalities=S.modalities; catch, warning('Using default Neuromag modalities'); S.modalities={'MEGPLANAR' 'MEGMAG'}; end;
try, S.pca_dim=S.pca_dim; catch, S.pca_dim=60; end; % rank to be use for the pca dimensionality reduction. If pca_dim=-1 then spm_pca_order is used to estimate pca_dim. If S.do_reduce_pca=1 then this is used prior to any source recon, if S.do_reduce_pca=0 then this is only used to determine the rank of the beamformer data covariance
try, S.regpc=S.regpc; catch, S.regpc=0; end; % data covariance regularisation
try, S.gridstep=S.gridstep; catch, S.gridstep=7; end; % space between dipoles in beamformer grid in mm
try, S.mask_fname=S.source_recon_mask_fname; catch, end;
try, S.useheadshape = S.useheadshape; catch, S.useheadshape=1; end
try, S.forward_meg = S.forward_meg; catch, S.forward_meg='MEG Local Spheres'; end % MEG head forward model set to 'Single Shell' or 'MEG Local Spheres'
try, S.mri = S.mri; catch, S.mri=''; end; % file path for MRI struct
try, S.fid_label = S.fid_label; catch, warning('Using default Neuromag fid labels'); S.fid_label.nasion='Nasion'; S.fid_label.lpa='LPA'; S.fid_label.rpa='RPA'; end; % To see what these should be look at: D.fiducials; fidnew.fid.label
try, S.fid_mnicoords = S.fid_mnicoords; catch, warning('Using default Neuromag fid std brain MNI coords');  end; % fids in mni coords, e.g. fid_mnicoords.nasian=[4.1 84.1 -11.7];fid_mnicoords.rpa=[83.8 -9.1 -65.1];fid_mnicoords.lpa=[-84.6 -3.3 -66.8]; Leave field empty to use standard brain fid coords
try, S.do_plots = S.do_plots; catch, S.do_plots=1; end;
try, S.work_in_pca_subspace = S.work_in_pca_subspace; catch, S.work_in_pca_subspace=1; end;
try, S.session_name = S.session_name; catch, S.session_name='session'; end;
try, S.force_pca_dim = S.force_pca_dim; catch, S.force_pca_dim=0; end;
try, S.reduce_rank = S.reduce_rank; catch, S.reduce_rank = 2; end;
try, S.type = S.type; catch, S.type = 'Scalar'; end;
try, S.hmm_av_class_occupancy = S.hmm_av_class_occupancy; catch, S.hmm_av_class_occupancy = 20; end;

try, report = S.report; catch, report = osl_report_setup([S.dirname '/report'],['Source recon (epoched)']); end;

if(isfield(S,'D')),
    error('Do not specify S.D. Specify S.D_continuous or S.D_epoched instead.');
end;

inverse_method=S.method;

source_recon_results=[];

clear matlabbatch;

conditions=S.conditions;

if strcmp(S.modalities{1},'EEG')
    modality_meeg='EEG';
else
    modality_meeg='MEG';
end

% check SPM MEEG objects passed in
only_epoched_data_provided=0;
if ~isempty(S.D_epoched) && isempty(S.D_continuous),
    only_epoched_data_provided=1;
    S.D=S.D_epoched;
elseif ~isempty(S.D_continuous),
    S.D=S.D_continuous;
else
    error('Need to specify S.D_continuous or S.D_epoched.');
end;

% do epoching if needed, but only if any epoch info is provided
do_epoching=0;
if ~only_epoched_data_provided
    if ~isempty(S.D_epoched) || ~isempty(S.epochinfo),
        do_epoching=1;
    end;
end;

if strcmp(S.forward_meg,'MEG Local Spheres') && S.work_in_pca_subspace,
    error('MEG Local Spheres not compatible with working in PCA subspace');
end;

if (strcmp(S.modalities{1},'MEGMAG') || strcmp(S.modalities{1},'MEGPLANAR')) && length(S.modalities)==1 && ~S.work_in_pca_subspace,
    error('Must work in PCA subspace if using only one MEG modality type, when there is more than one available.');
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%

disp(['Running beamformer for ' S.D]);
disp(['Will be designated ' S.session_name]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% move copy of SPM MEEG object to work on into the working directory
warning off;
mkdir(S.dirname);
warning on;

% move the spm file(s), giving it a unique name in the parent directory
S2=[];
S2.D=S.D;
tempstring = tempname;
tempstring = tempstring(end-12:end);
S2.newname=[S.session_name tempstring '_spm_meeg'];
[p spmname e] = fileparts(S2.newname);                               
S2.newname=[spmname '.mat'];   

S2.updatehistory=0;
D = spm_eeg_copy(S2);
runcmd(['mv ' D.path '/' D.fname ' ' S.dirname]);
runcmd(['mv ' D.path '/' D.fnamedat ' ' S.dirname]);
spm_filename=[S.dirname '/' D.fname];
D=spm_eeg_load(spm_filename);

% copy the file in the analysis directory in order to give it the standard
% name
S2=[];
S2.D = D;
S2.newname = [S.session_name '_spm_meeg'];
S2.updatehistory=0;
D = spm_eeg_copy(S2);

% delete the temporary files
matfn = fullfile(S.dirname,[S.session_name tempstring '_spm_meeg.mat']);
datfn = fullfile(S.dirname,[S.session_name tempstring '_spm_meeg.dat']);
runcmd(['rm ' matfn]);
runcmd(['rm ' datfn]);

spm_filename = fullfile(S.dirname,[S.session_name '_spm_meeg']);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% do bandpass filtering
if(length(S.freq_range)>1 && S.freq_range(2)>S.freq_range(1)),
    disp('Temporal filtering...');
    S2=[];
    S2.D=spm_filename;
    S2.use_fft_bandpass=1;
    
    if(0)
        
        % implement bandpass as a low pass followed by a high pass
        S2.filter.PHz=S.freq_range(2);
        S2.filter.band='low';
        Dnew = spm_eeg_filter_v2(S2);
        
        D.delete;
        D=Dnew;
        
        S2.D=D;
        S2.filter.PHz=S.freq_range(1);
        S2.filter.band='high';
        Dnew = spm_eeg_filter_v2(S2);
        
    else
        
        S2.filter.PHz=S.freq_range;
        S2.filter.band='bandpass';
        S2.use_fft_bandpass=1;
        Dnew = spm_eeg_filter_v2(S2);
        
    end;
    
    if S.bandstop_filter_mains
        % do notch filtering to remove mains noise
        notchBands = [48 52; 98 102; 148 152; 198 202; 248 252; 298 302; 348 352; 398 402; 448 452; 498 502];
        lowIn      = notchBands(:,1) > S.freq_range(1);
        highIn     = notchBands(:,2) < S.freq_range(2);
        doNotch    = find(lowIn & highIn);
        
        for iNotch = 1:numel(doNotch)
            S3              = [];
            S3.D            = Dnew;
            S3.filter.PHz   = notchBands(doNotch(iNotch),:);
            S3.filter.dir   = 'twopass';
            S3.filter.band  = 'stop';
            disp(['Band-stop filtering from ' num2str(S3.filter.PHz(1)) 'Hz to ' num2str(S3.filter.PHz(2)) 'Hz.']);
            Dnew = spm_eeg_filter_v2(S3);
        end % for iNotch = 1:numel(doNotch)
        
    end % if S.notch_filter_mains
    
    
    D.delete;
    D=Dnew;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% do epoching

% epochinfo will be in S.epochinfo, if its empty then look in D_epoched for D_epoched.epochinfo
if(do_epoching),
    disp('Epoching...');
    
    S2 = [];
    S2.D = D;
    
    if(~isempty(S.epochinfo)),
        S2.epochinfo=S.epochinfo;
    else
        try,
            D_epoched=spm_eeg_load(S.D_epoched);
            S2.epochinfo=D_epoched.epochinfo;
        catch
            error('No epoch info available');
        end;
    end;
    S2.epochinfo.padding = 0;
    S2.save=0;
    S2.reviewtrials=0;
    S2.bc=0;
    Dnew = spm_eeg_epochs(S2);
    
    D.delete;
    D=Dnew;
    
    %% get bad channels and trials from passed in S.D_epoched
    if(~isempty(S.D_epoched)),
        D_epoched_passed_in=spm_eeg_load(S.D_epoched);
        D = reject(D, 1:length(D.conditions), D_epoched.reject);
        if(length(D_epoched_passed_in.badchannels)>0),
            D = badchannels(D, D_epoched_passed_in.badchannels, ones(length(D_epoched_passed_in.badchannels),1));
        end;
        D.save;
    end;
    
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Establish time windows of interest

if isempty(S.time_range),
    S.time_range = [D.time(1) D.time(end)];
end
time_range=S.time_range;

Events = D.events; GoodEpochs=[D.time(1) D.time(end)];

if(D.ntrials==1),
    
    % if (ultimately) working on continuous data then use passed in time
    % range (which applies to continuous time)
    % BUT if (ultimately) working on continuous data, then work with full continuous data
    % and then apply passed in time range after epoching later on.
    if(do_epoching)
        time_range = [D.time(1) D.time(end)];
    end;
    
    % if continuous data then look for bad epochs and set time range accordingly
    
    for ev = 1:numel(Events)
        if isfield(Events,'type') && strcmp(Events(ev).type,'BadEpoch')
            GoodEpochs(end,2)=Events(ev).time;
            GoodEpochs(end+1,1)=Events(ev).time+Events(ev).duration;
            GoodEpochs(end,2)=D.time(end);
        end
    end
    
end

good_samples=zeros(1,D.nsamples);
for i=1:size(GoodEpochs,1),
    good_samples(D.indsample(GoodEpochs(i, 1)):D.indsample(GoodEpochs(i, 2)))=1;
end

samples_of_interest=zeros(1,D.nsamples);
for i=1:size(time_range,1),
    samples_of_interest(D.indsample(time_range(i, 1)):D.indsample(time_range(i, 2)))=1;
end

samples2use = samples_of_interest & good_samples;
woi=[D.time(find(diff([0 samples2use])==1))' D.time(find(diff([samples2use 0])==-1))'];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Establish trials

if strcmp(S.conditions{1},'all'),
    trials = D.pickconditions(D.condlist);
else
    
    trials = D.pickconditions(conditions);
    
    %     else
    %         trials=[];   % changed by DM
    %         for k=1:length(conditions)
    %             trials = [trials D.pickconditions(conditions{k})];
    %         end
    %     end
    %
    if isempty(trials)
        error('No trials matched the selection, check the specified condition labels');
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Establish setup channels

chanindmeg = strmatch(modality_meeg, D.chantype);

chanind = setdiff(chanindmeg, D.badchannels);
if isempty(chanind)
    error(['No good ' modality_meeg ' channels were found.']);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% save a copy of raw tra matrix before it is effected by spm_eeg_montage
% call when mapping onto pca space

if strcmp(modality_meeg,'MEG')
    D = save_raw_tra_to_D(D);
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% DO CO-REGISTRATION AND HEAD MODEL
% Need to do this before any montages (e.g. normalising), to ensure that
% neuromag_planar_baseline_corrections are applied correctly

disp('Co-registration and setting up forward model...');

S2=[];
S2.D = [D.path '/' D.fname];    % requires .mat extension
S2.mri=S.mri;
S2.useheadshape=S.useheadshape;
S2.forward_meg=S.forward_meg;
S2.fid_label=S.fid_label;
if(isfield(S,'fid_mnicoords')),
    S2.fid_mnicoords=S.fid_mnicoords;        
end;

S2.neuromag_planar_baseline_correction=S.neuromag_planar_baseline_correction;

% turned off display calls in:
% spm_cfg_eeg_inv_headmodel , spm_eeg_inv_datareg, spm_eeg_inv_datareg_ui,
% spm_eeg_inv_mesh_ui, spm_eeg_inv_icp

D=osl_forward_model(S2);

%% CHECK REGISTRATION AND PRODUCE DIAGNOSTIC PLOTS
if(1)
    
    fig_handles=[];
    fig_handles(1)= sfigure;
    
    set(fig_handles(1),'Position',[300 300 1200 800]);
    set(fig_handles(1),'Color', [1,1,1]);
    snugplot(1,3,1)
    spm_eeg_inv_checkdatareg_3Donly(D);
    view(-180,0)
    %title(['concatMefsession' num2str(counter) '_spm_meeg'])
    snugplot(1,3,2)
    spm_eeg_inv_checkdatareg_3Donly(D);
    view(-270,0)
    snugplot(1,3,3)
    spm_eeg_inv_checkdatareg_3Donly(D);
    view(130,18)

    report=osl_report_set_figs(report,'Coregistration',fig_handles);
    report=osl_report_print_figs(report);
           
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Normalise modalities using smallest eigenvalues
%% calculated using good channels and good trials, and over all woi
disp('Establish dimensionality and Normalising modalities...');

S2=S;
S2.D = D;
S2.samples2use=samples2use;
S2.trials=trials;
S2.do_plots=1;
S2.normalise_method=S.normalise_method;

%[ Dnew pcadims pcadim ] = normalise_sensor_data( S2 );

[ Dnew pcadims tmp tmp fig_handles fig_names] = normalise_sensor_data( S2 );

% set pcadim to min:
pcadim=min(pcadims);

% diagnostic plot of design matrix    
report=osl_report_set_figs(report,fig_names,fig_handles);
report=osl_report_print_figs(report);

D.delete;
D=Dnew;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% DO HMM

if(isfield(S,'hmm_num_states') && S.hmm_num_states~=0),
    
    %%
    if(isfield(S,'hmm_block')),
        block=S.hmm_block;
        
        if(S.hmm_num_states<0)
            NK=size(max(block.q_star),2);
        else
            NK=S.hmm_num_states;
            if size(max(block.q_star),2)~=NK,
                error('Incompatible number of states');
            end;
        end;
        
    else
        
        disp('Doing HMM...');
        
        tres=D.time(2)-D.time(1);
        tt=tres*length(find(samples2use))*length(trials);
        
        recommend_NK=min(S.hmm_pca_dim-15,round(tt/S.hmm_av_class_occupancy));
        disp(['There is ' num2str(tt) 'secs of data, recommend using no more than ' num2str(recommend_NK) ' HMM states.']);
        
        if(S.hmm_num_states<0)
            NK=recommend_NK;
        else
            NK=S.hmm_num_states;
        end;
        
        disp(['Inferring ' num2str(NK) ' HMM states']);
        
        Sh=[];
        
        Sh.data=reshape(D(chanind,find(samples2use),trials),length(chanind),(numel(find(samples2use))*length(trials)))';
        
        if(S.hmm_pca_dim>0)
            [allsvd,Apca]=pca(Sh.data,S.hmm_pca_dim);
            pinvApca=pinv(Apca);
            Sh.data=(pinvApca*Sh.data')';
        else
            error('hmm_pca_dim unspecified');
        end;
        
        Sh.data=normalise(Sh.data);
        
        Sh.NK=NK;
        
        Sh.num_starts=S.hmm_num_starts;
        
        [ hmm, block ]=run_multistart_hmm(Sh);
        
        NK=hmm.K;
        
    end;
    
    source_recon_results.block=block;
    
    hmm_class=zeros(1,size(D,2),size(D,3));
    hmm_class(1,find(samples2use),trials)=reshape(source_recon_results.block(1).q_star,[1,numel(find(samples2use)),length(trials)]);
    
    hmm_class_probs=zeros(NK,size(D,2),size(D,3));
    hmm_class_probs(:,find(samples2use),trials)=reshape(source_recon_results.block(1).gamma',[NK,numel(find(samples2use)),length(trials)]);
    
    %% plot results
    if(1),
          
        S2=[];
        S2.hmm=hmm;;
        S2.block=block;
        S2.tres=tres;
        S2.Apca=Apca;
        S2.NK=NK;
        
        [fig_handles fig_names fig_titles]=plot_hmm(S2);

        report=osl_report_set_figs(report,fig_names,fig_handles,fig_titles);        
        report=osl_report_print_figs(report);
 
        % if epoched data, plot epoched state courses
        if D.ntrials>1                    
        
            tmp=zeros(NK,size(hmm_class,2));
            leg={};
            for ii=1:NK,
                tmp(ii,:)=sum(hmm_class==ii,3);
                leg{ii}=num2str(ii);
            end;

            fig_handle=sfigure;
            plot(D.time(samples2use),tmp(:,samples2use)','LineWidth',2);
            plot4paper('time (s)','# of trials');
            legend(leg,'Location','NorthWest','FontSize',12); 
                         
            report=osl_report_set_figs(report,'hmm_epoched_state_occurences',fig_handle,'HMM Epoched State Occurrences');        
            report=osl_report_print_figs(report);
            
        end;
        
        if(0)
            data=normalise(Sh.data);
            sqdata=mean(abs(data),2);
            rng=max(sqdata)-min(sqdata);
            hmm_class2=hmm_class(1,find(samples2use),trials);
            hmm_class2=reshape(hmm_class2,size(hmm_class2,1),size(hmm_class2,2)*size(hmm_class2,3));

            ii=0;
            tmp=sqdata;
            fig_handle=sfigure;
            plot(tmp+rng*(ii) ,'b');
            ho;
            for ii=1:NK,

                tmp=sqdata;
                tmp(hmm_class2~=ii,:)=0;
                plot(tmp+rng*(ii) ,'b');

            end;
        end;
        
    end;
    
    clear Sh;
    
else,
    
    NK=1;
    
    hmm_class=zeros(1,size(D,2),size(D,3));
    hmm_class_probs=zeros(1,size(D,2),size(D,3));
    
    hmm_class(1,find(samples2use),trials)=1;
    hmm_class_probs(1,find(samples2use),trials)=1;
end;


% make a new state for eyeblinks/artefacts, and remove it from the other states
if isfield(S, 'artefact_chanlabel') && ~isempty(S.artefact_chanlabel)
    artefactchanind = find(strcmp(S.artefact_chanlabel,D.chanlabels));
    if isempty(artefactchanind)
        error('Can''t find the specified artefact channel in the D-object!');
    end
    maxclass = max(unique(hmm_class));
    artfctclasssinds = logical(D(artefactchanind,:,:));
    hmm_class(artfctclasssinds) = maxclass + 1;
    % the highest indexed state is now the artefact state
end

% add channel
Sc=[];
Sc.D=D;
Sc.newchandata=[hmm_class; hmm_class_probs];
Sc.newchanlabels{1}='Class';
Sc.newchantype{1}='CLASS';

for kk=1:NK,
    Sc.newchanlabels{kk+1}=['ClassPr' num2str(kk)];
    Sc.newchantype{kk+1}='CLASSPR';
end;

[ Dnew ] = osl_concat_spm_eeg_chans( Sc );

D.delete;
D=Dnew;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% do dimensionality reduction

% use montage to move everything into PCA subspace from now on
if(S.work_in_pca_subspace)
    disp('Doing PCA dimensionality reduction...');
    
    if strcmp(modality_meeg,'EEG'),
        error('Not currently implemented');
    end;
    
    if(S.regpc~=0)
        warning('No covariance regularisation is done if working in PCA subspace.');
    end;
    
    S2=[];
    
    S2.D=D;
    S2.modality_meeg=modality_meeg;
    S2.method='pca_adapt';
    
    S2.settings.force_pca_dim=1;
    S2.settings.ncomp=pcadim;
    S2.conditions=conditions;
    S2.woi = woi;
    S2.channels=chanind;
    
    [Dnew pca_montage] = spm_eeg_reduce(S2);
    
    source_recon_results.pca_tra=pca_montage.tra;
    
    D.delete;
    D=Dnew;
    
    % update the modality names
    if strcmp(modality_meeg,'MEG')
        S.modalities={'MEGPCACOMP'};
        pcadim = sum(strcmp(D.chantype,'MEGPCACOMP'));
    else
        S.modalities={'EEGPCACOMP'};
        pcadim = sum(strcmp(D.chantype,'EEGPCACOMP'));
    end;
    
end;

clear dat;

disp(['PCA rank from now on is: ' num2str(pcadim)]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% DO CO-REGISTRATION AND HEAD MODEL
% Need to at least re-do the head model following doing the montage

disp('Co-registration and setting up forward model...');

S2=[];
S2.D = [D.path '/' D.fname];    % requires .mat extension
S2.mri=S.mri;
S2.useheadshape=S.useheadshape;
S2.forward_meg=S.forward_meg;
S2.fid_label=S.fid_label;
S2.neuromag_planar_baseline_correction='none';

% turned off display calls in:
% spm_cfg_eeg_inv_headmodel , spm_eeg_inv_datareg, spm_eeg_inv_datareg_ui,
% spm_eeg_inv_mesh_ui, spm_eeg_inv_icp

D=osl_forward_model(S2);

if(S.do_plots)
    %spm_eeg_inv_checkmeshes(D);
    spm_eeg_inv_checkdatareg(D);
    %spm_eeg_inv_checkforward(D, 1);
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% do beamformer batch stuff

disp('Beamforming...');

disp(['Source localization using ' inverse_method]);  % changed by DM

i=0;

%%%%%%%%%%%%%
%%
i=i+1;

matlabbatch{i}.spm.tools.beamforming.data.D{1} = [D.path '/' D.fname];
matlabbatch{i}.spm.tools.beamforming.data.dir = {S.dirname};
matlabbatch{i}.spm.tools.beamforming.data.val = 1;
if strcmp(modality_meeg,'MEG')  % added by DM
    matlabbatch{i}.spm.tools.beamforming.data.space = 'Head';
else
    matlabbatch{i}.spm.tools.beamforming.data.space = 'Native';
end

%%%%%%%%%%%%%
%%
i=i+1;

%matlabbatch{i}.spm.tools.beamforming.sources.BF(1) = cfg_dep('Prepare data: BF.mat file', substruct('.','val', '{}',{i-1}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','BF'));

matlabbatch{i}.spm.tools.beamforming.sources.BF(1) = cfg_dep;
matlabbatch{i}.spm.tools.beamforming.sources.BF(1).tname = 'BF.mat file';
matlabbatch{i}.spm.tools.beamforming.sources.BF(1).tgt_spec{1}(1).name = 'filter';
matlabbatch{i}.spm.tools.beamforming.sources.BF(1).tgt_spec{1}(1).value = 'mat';
matlabbatch{i}.spm.tools.beamforming.sources.BF(1).tgt_spec{1}(2).name = 'strtype';
matlabbatch{i}.spm.tools.beamforming.sources.BF(1).tgt_spec{1}(2).value = 'e';
matlabbatch{i}.spm.tools.beamforming.sources.BF(1).sname = 'Prepare data: BF.mat file';
matlabbatch{i}.spm.tools.beamforming.sources.BF(1).src_exbranch = substruct('.','val', '{}',{i-1}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1});
matlabbatch{i}.spm.tools.beamforming.sources.BF(1).src_output = substruct('.','BF');

matlabbatch{i}.spm.tools.beamforming.sources.plugin.mni_coords.pos = S.mni_coord;
matlabbatch{i}.spm.tools.beamforming.sources.reduce_rank = [S.reduce_rank 3];

%%%%%%%%%%%%%
%%

i=i+1;
%matlabbatch{i}.spm.tools.beamforming.features.BF(1) = cfg_dep('Define sources: BF.mat file', substruct('.','val', '{}',{i-1}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','BF'));

matlabbatch{i}.spm.tools.beamforming.features.BF(1) = cfg_dep;
matlabbatch{i}.spm.tools.beamforming.features.BF(1).tname = 'BF.mat file';
matlabbatch{i}.spm.tools.beamforming.features.BF(1).tgt_spec{1}.name = 'filter';
matlabbatch{i}.spm.tools.beamforming.features.BF(1).tgt_spec{1}.value = 'mat';
matlabbatch{i}.spm.tools.beamforming.features.BF(1).sname = 'Define sources: BF.mat file';
matlabbatch{i}.spm.tools.beamforming.features.BF(1).src_exbranch = substruct('.','val', '{}',{i-1}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1});
matlabbatch{i}.spm.tools.beamforming.features.BF(1).src_output = substruct('.','BF');

matlabbatch{i}.spm.tools.beamforming.features.whatconditions.condlabel=conditions;
for jj=1:length(S.modalities),
    matlabbatch{i}.spm.tools.beamforming.features.whatmodalities.modalitylabel{jj}=S.modalities{jj};
end;

matlabbatch{i}.spm.tools.beamforming.features.woi = woi;
matlabbatch{i}.spm.tools.beamforming.features.plugin.regmulticov.lambda = S.regpc;
if(S.regpc~=0)
    error('S.regpc~=0 currently not supported');
end;

%%%%%%%%%%%%%
%%

i=i+1;
%matlabbatch{i}.spm.tools.beamforming.inverse.BF(1) = cfg_dep('Covariance features: BF.mat file', substruct('.','val', '{}',{i-1}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','BF'));

matlabbatch{i}.spm.tools.beamforming.inverse.BF(1) = cfg_dep;
matlabbatch{i}.spm.tools.beamforming.inverse.BF(1).tname = 'BF.mat file';
matlabbatch{i}.spm.tools.beamforming.inverse.BF(1).tgt_spec{1}.name = 'filter';
matlabbatch{i}.spm.tools.beamforming.inverse.BF(1).tgt_spec{1}.value = 'mat';
matlabbatch{i}.spm.tools.beamforming.inverse.BF(1).sname = 'Covariance features: BF.mat file';
matlabbatch{i}.spm.tools.beamforming.inverse.BF(1).src_exbranch = substruct('.','val', '{}',{i-1}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1});
matlabbatch{i}.spm.tools.beamforming.inverse.BF(1).src_output = substruct('.','BF');

if strcmp(inverse_method,'beamform')
    matlabbatch{i}.spm.tools.beamforming.inverse.plugin.lcmv_multicov.pca_order = pcadim;
    matlabbatch{i}.spm.tools.beamforming.inverse.plugin.lcmv_multicov.type = S.type;
elseif strcmp(inverse_method,'mne_diag_datacov')
    mne_lambda=1;
    matlabbatch{i}.spm.tools.beamforming.inverse.plugin.mne_multicov.noise_cov_type='diag_datacov';
    matlabbatch{i}.spm.tools.beamforming.inverse.plugin.mne_multicov.lambda = mne_lambda;
    matlabbatch{i}.spm.tools.beamforming.inverse.plugin.mne_multicov.type = S.type;
elseif strcmp(inverse_method,'mne_datacov')
    mne_lambda=1;
    matlabbatch{i}.spm.tools.beamforming.inverse.plugin.mne_multicov.noise_cov_type='datacov';
    matlabbatch{i}.spm.tools.beamforming.inverse.plugin.mne_multicov.lambda = mne_lambda;
    matlabbatch{i}.spm.tools.beamforming.inverse.plugin.mne_multicov.type = S.type;
elseif strcmp(inverse_method,'mne_eye')
    mne_lambda=1;
    matlabbatch{i}.spm.tools.beamforming.inverse.plugin.mne_multicov.noise_cov_type='eye';    
    matlabbatch{i}.spm.tools.beamforming.inverse.plugin.mne_multicov.lambda = mne_lambda;
    matlabbatch{i}.spm.tools.beamforming.inverse.plugin.mne_multicov.type = S.type;
elseif strcmp(inverse_method,'champagne')
    matlabbatch{i}.spm.tools.beamforming.inverse.plugin.champagne_multicov.type = S.type;
else
    disp('inversion method unknown!');
end

%% call jobman

spm_jobman('run', matlabbatch);
source_recon_results.BF=load([S.dirname '/BF.mat']);

source_recon_results.pca_order = pcadim;
source_recon_results.type = S.type;

source_recon_results.woi=woi;
source_recon_results.samples2use=samples2use;
source_recon_results.trials=trials;
source_recon_results.chanind=chanind;
source_recon_results.report=report;

runcmd(['rm -rf ' S.dirname '/BF.mat']);

end